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A  mathematical  simulation  model  was  developed  that  can  determine  the  three-dimensional  wind  veloc¬ 
ity  field  over  a  complex  terrain.  The  Tenes  area  in  the  Valley  of  Cheliff  in  Algeria  was  used  as  a  case  study. 
This  region  is  exposed  to  south-west  circulation  that  makes  it  favorable  to  the  use  of  wind  energy.  Knowl¬ 
edge  of  wind  fields  is  crucial  for  predicting  the  dispersion  of  pollutants,  for  forecasting  meteorological 
weather,  for  fire  spread  prediction  and  in  the  design  and  implementation  of  wind  turbines.  By  means 
of  a  mass  consistent  model,  an  in-house  program  was  developed  to  calculate  the  three-dimensional 
wind  velocity  field  in  the  study  region.  The  model  was  supported  by  a  numerical  box  in  which  flow 
through  is  allowed  for  in  the  upper  and  lateral  boundaries.  The  bottom  boundary  through  which  no  flow 
through  occurs  was  determined  by  the  topographic  relief  at  the  surface.  From  measured  wind  veloci¬ 
ties,  observed  values  were  calculated  by  interpolation-extrapolation.  Using  an  optimization  method,  the 
adjusted  velocities  were  obtained  from  constraints,  observed  velocities  and  the  continuity  equation.  The 
model  was  verified  with  wind  point  data,  the  relative  error  did  not  exceed  6%. 
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1.  Introduction 
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The  interest  in  wind  energy  production  has  increased  consider¬ 
ably  due  to  improvements  in  turbine  technology  and  the  growing 
need  for  power  production  from  renewable  sources  [1-3].  One  of 
the  most  important  issues  to  deal  with  is  the  choice  of  an  adequate 
site  for  a  wind  farm  [4].  Different  methodologies  have  been  devel¬ 
oped  to  investigate  and  define  the  wind  climate  of  regions  using 
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Nomenclature 

X  Lagrange  multiplier 

variances  of  the  observed  field 
ot[  Gauss  precision  module 

F  objective  function 

V  Computational  domain 

u°,v°  and  w°  fixed  observed  velocity  value  in  x,y  and  z  direc¬ 
tion 

u ,  v  and  w  velocity  adjusted  value  in  x ,  y  and  z  direction 
Ax,  Ay,  Az  step  in  x,  y  and  z  direction 


actual  measured  data  [5].  Mortensen  and  Petersen  [6]  for  exam¬ 
ple  have  shown  that  numerical  models  are  suitable  for  simulating 
wind  fields.  Such  models  can  indicate  the  best  sites  for  installation 
of  wind  turbines. 

Due  to  a  boost  in  the  world’s  energy  consumption  as  a  result 
of  population  growth  and  industrialization,  renewable  energy  (i.e. 
solar,  wind,  geothermal  and  wave)  has  increased  in  importance  [7]. 
A  study  by  Ludwig  and  Byrd  [8]  has  shown  that  meteorological  and 
air  pollution  problems  frequently  require  the  interpolation  of  com¬ 
plete,  three-dimensional  wind  field  flows  from  a  limited  number  of 
observations  at  specific  location.  Modeling  atmospheric  flow  is  dif¬ 
ficult,  especially  when  simulating  wind  fields  over  irregular  terrain 
[9], 

Different  types  of  wind  field  models  have  been  developed  in 
the  last  few  years  [10].  They  can  be  classified  into  two  main  mod¬ 
els:  diagnostic  (kinetic)  and  prognostic  (dynamic).  Kinetic  models 
are  those  capable  of  reconstructing  a  steady-state  wind  field  start¬ 
ing  from  a  set  of  initial  experimental  data.  Their  equations  contain 
no  time  dependent  terms.  On  the  other  hand,  dynamic  models 
describe  the  time  evolution  of  meteorologically  variable  fields 
starting  from  an  initial  state  and  merging  the  effects  of  possible 
time  variations  of  boundary  conditions. 

Diagnostic  (i.e.  kinetic)  models  follow  two  main  approaches 
[11]:  a  linearized  model  which  is  a  simplified  steady-state  solu¬ 
tion  of  the  Navier-Stokes  equation,  and  the  mass-consistent  model 
which  analyses  available  meteorological  data,  with  physical  con¬ 
straints  like  mass-conservation.  It  is  a  general  practice  to  use  mass 
consistent  models  for  simulation  of  stationary  three-dimensional 
wind  fields  in  complex  terrain.  The  advantage  of  a  mass  consistent 
model  compared  with  primitive  equation  models  is  the  relatively 
short  computing  time.  However,  as  a  result  it  is  not  possible  to  look 
at  complex  phenomenon  such  as  turbulence  [12]. 

The  implementation  of  a  mass  consistent  wind  field  is  typically 
undertaken  using  available  meteorological  data.  The  procedure  has 
been  effective  for  predicting  winds  for  emergency  response  con¬ 
taminant  dispersion,  for  weather  forecasts,  in  wind  energy  studies 
and  for  storm  flooding  prediction  [13].  Furthermore,  some  applica¬ 
tions  require  that  the  derived  wind  field  be  non  divergent. 

The  mass-consistent  models  reconstruct  3-D  wind  fields 
starting  from  vertical  profiles  and  near-ground  wind  measure¬ 
ments  through  a  two-step  procedure  [14].  The  observed  data 
needed  for  the  mass  consistent  model  are  provided  by  an 
interpolation-extrapolation  scheme  using  information  available  at 
a  given  site  to  determine  velocity  components  at  each  grid  point 
above  the  topography.  Typically  the  available  simultaneous  data 
consists  of  several  horizontal  wind  speed  and  direction  measure¬ 
ments  near  the  earth’s  surface,  a  vertical  profile  of  horizontal  wind 
speed  and  direction  and  synoptic  analysis. 

The  wind  data  interpolation  method  plays  a  crucial  role  in  deter¬ 
mining  the  final  wind  field  features.  Different  wind  fields  can  be 
generated  starting  from  the  same  data  and  changing  only  the  inter¬ 
polation  method.  The  adjustment  procedure  is  generally  obtained 


Table  1 

Characteristics  of  the  measuring  station  [19]. 


Site 

Longitude  (°) 

Latitude  (°) 

Altitude  (M) 

Period  of 

measure 

(year) 

Height  of 
mast  (m) 

Tenes 

1.33 

36.5 

17 

5 

10 

by  employing  the  variational  analysis  technique  [15].  Early  research 
was  conducted  by  Sasaki  [16-18]  who  developed  the  first  theo¬ 
retical  model  for  the  calculation  of  wind  speed.  Many  different 
mass-consistent  models  have  been  developed  starting  from  the 
late  seventies  [19-29].  Various  other  works  relating  the  results  of 
air  pollution  dispersion  modeling  in  complex  terrain  are  entirely 
dependent  on  adequate  descriptions  of  the  three-dimensional  wind 
field  [30].  Furthermore,  as  part  of  modernization  efforts  of  the 
Atmospheric  Release  Advisory  Capability  (ARAC)  project,  Chan  and 
Sugiyama  [31]  developed  a  new  diagnostic  model  for  generating 
mass  consistent  wind  fields  over  a  continuous  terrain.  In  addition, 
Leone  et  al.  [32]  have  used  wind  fields  to  drive  ARAC’s  new  disper¬ 
sion  model.  This  model  is  going  to  replace  the  current  operational 
code  described  by  Sherman  [19]  which  employs  stair-step  topogra¬ 
phy  and  constant  grid  spacing.  The  topography  strongly  influences 
the  boundary  layer  flow  in  such  a  manner  that  the  geographic  vari¬ 
ations  either  can  modify  on  in  some  cases  generate  circulations  in 
conjunction  with  diurnal  heating  cycles  [33,34]. 

The  theoretical  basis  for  the  mass  consistent  model  was  devel¬ 
oped  by  Sasaki  [16-18].  The  general  variational  analysis  formalism 
defines  an  integral  function  whose  external  solution  minimizes  the 
variance  of  the  difference  between  the  observed  and  analyzed  vari¬ 
able  values  subject  to  physical  constraints  which  are  the  continuity 
equation  and  the  measured  velocities. 

The  aim  of  this  paper  is  to  develop  a  mathematical  simulation 
model  that  can  determine  the  three-dimensional  wind  velocity 
field  over  a  complex  terrain  using  the  Valley  of  Cheliff  in  the  Tenes 
area  in  Algeria  as  a  case  study.  This  region  has  a  complex  relief  with 
a  set  of  mountains  and  valleys  with  strong  slopes.  In  the  adjustment 
model  for  wind  flow,  the  relief  is  analyzed  explicitly  with  the  use  of 
velocity  measurements  in  the  field  using  a  large  number  of  points 
in  a  three-dimensional  grid.  The  data  are  interpolated  and  extrap¬ 
olated  by  a  least  squares  method  on  the  grid  of  an  Euler  model.  The 
numerical  results  are  compared  with  the  measurement  values  at 
selected  points  on  the  real  terrain. 

2.  Methodology 

2.2.  Case  study  area 

The  study  area  is  located  in  the  town  of  Tenes  in  Algeria.  It  is  a 
coastal  city  which  is  partially  surrounded  by  the  Dahra  mountain 
range.  The  topography  of  this  region  is  complex  with  a  varia¬ 
tion  in  altitude  between  0  and  700m  (Fig.  la).  The  correspondent 
mathematical  topographic  model  is  given  in  Fig.  lb.  The  data  file 
that  defines  the  terrain  relief  by  the  three  coordinates  contains 
374,037  points.  Digital  processing  was  used  to  order  the  relief  struc¬ 
ture.  The  site  is  characterized  by  an  important  wind  potential  area 
[35].  The  weather  station  at  this  location  has  5  years  of  continu¬ 
ous  measurements  of  three-hourly  speed  and  wind  direction  [36]. 
Characteristics  of  the  measuring  station  are  given  in  Table  1. 


2.2.  Mathematical  model  formulation 

The  study  is  presented  as  an  optimization  problem  with  con¬ 
straint  conditions.  The  objective  function  is  given  by  an  integral 
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Fig.  1.  Digital  and  mathematical  relief  of  Tenes. 


functional  [16-18]: 


E=  /( F(u,v,w,  X)dxdydz 


(1) 


V 


where  F  is  a  function  expressed  by 

F(u,  v,  w,  X)  =  a^{u  -  u°f  +  a2(v  -  v°f  +  aj(w 

,  .  du  dv  dw 
+  X  (  —  +  —  +  — 
dx  dy  dz 


-  w0/ 


(2) 


where  u ,  v,  and  w  are  the  argument  functions  which  represent 
respectively  the  components  wind  velocity  in  x,y,  and  z  direction  of 
the  modeling  domain  V.  The  mass  consistent  model  constructs  an 
observed  initial  guess  of  the  wind  velocities  component  u°,  v°  and 
w°  over  the  entire  model.  The  initial  estimates  are  formed  from 
an  interpolation-extrapolation  scheme  that  uses  existing  surface 
wind  observations  contained  within  the  modeled  area.  The  values  u , 
v  and  w  are  obtained  by  adjusted  scheme  from  the  observed  veloc¬ 
ities  by  minimizing  the  functional  E.  The  parameter  X(x,  y ,  z)  is  the 
Lagrange  multiplier;  and  values  of  oti  are  Gauss  precision  moduli 
taken  to  be  a2  =  1  /2a2.  a*  is  an  observation  error  and/or  variances 
of  the  observed  field  from  the  desired  adjusted  field.  Identical  Gauss 
precision  moduli  are  assumed  for  the  horizontal  directions.  The 
ratio  a  =  (aq \cl2 )2  =  (a2 \o\  )2  is  related  to  stability  and  turbulence  of 
atmosphere  [19]. 

The  variational  calculus  gives  the  Euler-Lagrange  equation  for 
u ,  v  and  w,  respectively  [37]: 


where  n  is  the  outward  positive  unit  normal  to  the  boundary  of  vol¬ 
ume  V.  Eq.  (5)  implies  that  on  the  boundary,  either  of  the  following 
is  satisfied: 

Vu  =  0  or  A  =  0  (6) 


Eq.  (5)  represents  appropriate  conditions  for  “flow-through”  lateral 
and  upper  boundaries.  Either  the  normal  velocity  component  vari¬ 
ation  or  X  must  be  zero  on  a  boundary.  The  appropriate  condition 
for  “flow-through”  boundary  (upper  and  lateral  boundary)  is  X  =  0. 
In  this  case  the  derivate  of  X  in  the  direction  across  the  boundary  is 
in  general  not  zero.  The  lower  boundary  coincides  with  the  ground 
and  represents  a  “no-flow-through”  solid  boundary.  The  observed 
velocity  at  that  boundary  is  zero  [19,25,42]. 

Substituting  for  F by  putting  Eq.  (2)  into  Eq.  (3),  we  obtain: 


w  =  w 


1  dX 

+  2a2  dx 
1  dX 


+ 

o 


2  ol\  3y 

1  dX 

+  2  ol\  dz 


The  continuity  equation  is  given  by: 


du 

dx 


dv 


dw 

dz 


=  0 


The  equation  for  X  is  derived  by  differentiating  Eq.  (7)  and  sub¬ 
stituting  the  results  into  Eq.  (8)  giving: 


dF 

d  1 

f  dF 

du 

dx  ' 

y  dux 

dF 

d  1 

(  dF 

dv 

dy  ' 

K.dVx. 

dF 

d  1 

(  dF 

dw 

dz  ' 

i  dwx 

where 

=  0 

=  0 

=  0 


^  +  =  9^+9^  ,q) 

dx2  dy 2  \a2)  dz2  1  (  dx  dy  dz  )  1  ’ 

Eq.  (9)  is  an  elliptical  differential  equation  (Poisson’s  equation). 
It  will  be  solved  for  X  with  the  specific  boundary  conditions  of  Eq. 
(5).  The  adjusted  wind  velocity  components  are  calculated  using 
Eq.(7). 

2.3.  Method  of  resolution 


ux 


du 

dx 


dv 

dy 


ux 


dw 

dz 


Eq.  (3)  are  subjected  to  the  boundary  conditions: 


TD,yu =o 


This  case  study  can  be  described  as  a  constraint  optimized  prob¬ 
lem  where  Eq.  (1)  is  the  objective  function  and  the  constraint 
conditions  are  the  continuity  equation  and  the  measured  values 

(4)  in  the  volume  V.  Mass-consistent  optimization  is  used  to  solve  the 
problem.  The  basic  procedure  to  get  the  wind  field  velocities  by  the 
mass-consistent  method  is  given  by  the  following  steps: 

The  Lagrange  parameter  X^  is  determined  by  Eq.  (9)  by  the 
finite  difference  method  in  the  discretised  volume  V  which  is  a  rect¬ 
angular  box  set  on  the  topographic  relief.  The  grid  model  is  built 

(5)  with  a  choice  of  Ax  and  Ay  such  that  the  grid  takes  in  account 
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the  more  interesting  topographic  detail  features  of  the  site  (e.g. 
mountains,  canyons  and  valley  areas).  The  computational  volume 
is  shown  in  Fig.  1. 

By  using  a  central  difference  scheme  for  second  derivatives  for 
a  grid  point  (i j,/<),  Eq.  (9)  is  rewritten  as: 


^Mj,k  2kijfk  +  hj-l,j,k  +  ^ ij+l,k  2A ijk  +  X i,j-l,k 


+ 


Ax2 
1  \ 


\  X\  j  b  1 1  2Xi  j  b  +  X 


i,j,k  ~r  Ai,j,k-1 


at 


A  z2 


Ay2 


=  -20'2£q 


(10) 


With  the  divergence  £0  of  the  interpolated  wind  field  is  given  by: 
du°  dv°  9w° 


e°-lF  +  ^+  3 z 


(11a) 


Fig.  2.  Finite  difference  grid  of  the  constraint  problem. 

3.  Results  and  discussion 


Using  central  difference  for  the  first  derivate  £q  is  written  as: 


u°,  1  . ,  -  u9  1  w? . ,  , ,  -  w° . ,  1 

i+l ,j,k  i-l ,j,k  ij+l,k  i,J-l,k  .  i,j,k+ 1  i,J,k- 1 

F0  -  - ^ - 1 - - F 


2  Ax 


2  Ay 


2Az 


(lib) 


where  the  observed  values  u9 .  ,  v° . ,  and  w° . ,  were  determined  by 
least  square  weighted  method  from  the  measured  wind  speed.  The 
method  uses  the  three  nearest  measured  points  to  determine  the 
observed  wind  speed  at  each  grid  point  (ij,/<).  Once  the  observed 
velocities  are  determined,  £o  is  introduced  in  Eq.  (10).  A  set  of  dif¬ 
ference  equation  in  X{^k  is  now  obtained  at  each  interior  grid  point. 
The  entire  system  is  solved  simultaneously  for  the  value  of  X{^k.  The 
equation  set  is  solved  iteratively,  using  a  successive  over-relaxation 
method. 

In  altitude  z  the  observed  values  (either  measured  or  adjusted 
values)  are  extrapolated  with  the  following  equation  [38]: 

(,2) 


where  the  exponent  p  is  determined  by  atmospheric  stability  con¬ 
ditions  or  from  a  least  squares  fit  to  multiple  level  tower  data. 

Once  the  Lagrange  parameters  are  determined,  the  adjusted 
wind  field  values  are  calculated  using  central  difference  for  Xtj  k 
and  mean  value  for  observed  wind  value  of  Eq.  (7): 


l/fi  'in  n  \  1  ( A+lJ,k  A— lj,k 

w  j  k  =  -  (uu  1  . ,  +  2uu . .  +  uu  1  )  H - -  [  - — — - — 

lJ’K  4V  i j,k  i-\,j,kJ  2cl2x  y  2 Ax 

L  n  ~  n  n  s  1  /  ^-i,j+\,k  ~  ^i,j-\,k 

vu,k  =  4  (vij+i,k  +  2vi,j,k  +  yij-u)  +  ^2  y  2A y 

1  f  1  —  1 


(13) 


'•W  =  3K-.t+,  +2wy,k  +<„-,)+  ( 


2Az 


3.2.  Numerical  development  of  the  terrain  and  grid  generation 

By  letting  x,  y  represent  the  horizontal  direction  and  z  the  verti¬ 
cal  direction  over  a  modeling  domain  V,  the  data  are  arranged  in  the 
y  direction  and  subdivided  in  columns.  The  choice  of  the  grid  size 
is  often  determined  by  the  necessity  to  show  the  real  topographic 
features.  Fig.  2  shows  a  cross  section  of  the  computational  volume. 
The  topography  may  protrude  through  the  top  of  the  grid  poten¬ 
tially  producing  unconnected  regions  within  the  grid  volume.  This 
feature  can  be  used  to  study  wind  fields  in  valleys,  canyons  and 
mountain  areas.  The  numerical  relief  shown  in  Fig.  3  indicates  that 
the  area  of  interest  for  this  model  is  a  rectangular  box  set  on  the 
earth’s  surface  with  the  bottom  of  the  box  located  at  the  lowest 
topographic  point.  The  dimensions  of  the  box  are  determined  by 
the  application  requirements  and  computer  storage  limitations. 

The  computational  domain  is  meshed  by  18x6x10  nodes  with 
step  Ax  =  Ay  =  2000  m  and  Az  =  1 00  m,  as  shown  in  Fig.  3.  The  val¬ 
ues  of  o\  and  <r2  were  taken  to  be  1  and  0.01  ms-1,  respectively. 
The  wind  speed  and  direction  at  the  upper  boundary  of  the  grid 
were  temporally  interpolated  from  the  two  vertical  profile  obser¬ 
vations  and  assumed  to  be  spatially  constant.  The  horizontal  wind 
measurements  were  adjusted  to  a  height  of  50  m  above  the  ter¬ 
rain  using  a  power  law  exponent  of  0.2  [39].  The  constant  height 
interpolated  and  extrapolated  wind  field  is  processed  to  define  hor¬ 
izontal  wind  components  at  each  grid  point  about  the  topography 
boundary.  The  foregoing  applications  differ  in  terrain,  meteorolog¬ 
ical  data  availability,  meteorological  condition,  grid  resolution,  and 
convergence  criterion. 

3.2.  Observed  win  speed  and  Lagrange  parameters 


The  algorithm  for  calculating  the  adjusted  wind  speed  values  is 

determined  as  follows: 

•  Obtain  three-dimensional  topographic  data  of  the  complex  ter¬ 
rain. 

•  Develop  three-dimensional  discrete  reliefs  for  finite  difference 
method. 

•  Determine  3-D  numerical  box  for  finite  difference  method  appli¬ 
cation. 

•  Localize  measured  values  of  wind  velocities  on  the  discrete  bot¬ 
tom  boundary  of  the  box. 

•  Determine  the  observed  wind  velocities  (u°,  v°  and  w°)  by  least- 
squares  method. 

•  Calculate  the  Lagrange  parameters  in  the  volume  box  by  finite 
difference  method. 

•  Use  SOR  method  to  resolve  set  of  equations  compute 
adjusted  velocity  components  (u,  v  and  w)  within  the  box  by 
Euler-Lagrange  method  (solution  of  the  problem  by  variational 
form). 


The  measured  values  are  extracted  from  results  of  three- 
dimensional  observed  velocities  ( u9  ,  v9. , ,  w?. , )  calculated  by 
the  AILOS  software  [40].  Fig.  4  shows  11  measured  points  (black 
symbols)  which  were  chosen  on  the  site  at  10  m  from  the  surface. 

The  discretised  Eq.  (9)  is  solved  to  give  the  Lagrange  parameters 
values  X[j  k  at  each  node  of  the  3D-grid.  The  matrix  coefficients  of 
the  equation  system  are  diagonally  dominant  and  asymmetric.  The 
iterative  successive  over  relaxation,  SOR,  method  [37]  was  used  to 
solve  the  correspondent  equation  system.  Eq.  (9)  can  be  written  as: 


(1  -  co)  X 


it- 1 
i,j,k 


+  CO  • 


A  /  k  +  L  ,  ■  .  +  L.1 .  +  klt7\  +  +  lambda.  7\  +  B(Xlt7 }  +  klt7}  ) 

i+l ,j,k  i-lj,fc  ij-l,k  ij+l,k  v  i,J,k+ 1  ij,k-ly 


C 


(14) 


where 

Aij,k  =  2(“i  Ax)2e0 


(15) 
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Fig.  3.  Three-dimensional  grid  computational  volume. 
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Fig.  4.  Observed  velocities  (colored  lines)  calculated  by  AILOS  software  [40]  and  the 
measured  points  (11  black  symbols).  (For  interpretation  of  the  references  to  color 
in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  the  article.) 


difference  of  over  500  m  between  the  highest  and  the  lowest  layer 
was  observed.  The  model  complex  relief  is  composed  of  a  set  of 
mountains  or  hills  with  a  high  slope.  The  major  characteristics  of 
the  flow  over  a  hill  were  determined  not  only  by  the  hill  shape  but 
also  by  its  size  and  by  atmospheric  stability.  The  horizontal  wind 
measurements  shown  by  the  colored  contour  lines  in  Fig.  4  were 
taken  from  the  observed  velocities  calculated  by  AILOS  software  in 
the  study  region  of  Tenes.  The  accuracy  of  the  wind  field  is  depen¬ 
dent  on  the  specification  of  a.\  and  a2.  For  a  stable  atmosphere  the 
horizontal  adjustment  were  given  by  =  1  and  a2  =  10;  then  the 
adjustment  is  predominant  in  the  horizontal  component  instead  of 
vertical  component  since: 


(19) 


The  developed  computer  code,  applied  on  the  numerical  relief 
in  shown  in  Fig.  3,  produces  nondivergent  three-dimensional  wind 
at  1080  grid  points  with  the  use  of  a  Pentium  4,  CPU  2.66  GHz. 


B  =  a 


C  =  4 


1  +  a 


(16) 

(17) 


where  Ay  =  Ax,  the  coefficients  Aijtk,  B  and  C  have  known  values. 
The  variable  specifies  the  current  value  of  the  iteration  in  the  SOR 
method. 

The  relaxation  factor  co  was  determined  experimentally  and  was 
equal  to  1.78.  This  value  is  chosen  so  as  to  obtain  the  minimum 
number  of  iterations.  The  iterative  solution  is  considered  converged 
when  the  relative  change  of  every  A.^  is  less  than  a  prescribed 
value.  The  iterative  process  is  stopped  when  the  condition  on  the 
solution  accuracy  is  achieved 


A*  -  A*-1 
A J 


■<  £* 


(18) 


where  s*  is  a  given  accuracy  value. 


3.3. 1 .  Wind  field  velocity  in  horizontal  planes 

Fig.  5a  shows  the  horizontal  wind  field  velocity  at  plane  (x,  y) 
above  the  topographic  relief.  The  wind  values  were  adjusted  with 
a  step  height  of  100  m  in  the  z  direction  using  Eq.  (12).  The  expo¬ 
nent  p  has  a  value  of  0.2  in  a  period  where  the  atmosphere  was 
neutrally  steady.  The  highest  wind  velocities  were  observed  in  the 
intervals  [364,366]  [352,356]  in  the  x  direction.  An  interesting  area 
is  a  valley  with  a  higher  depth  from  the  reference  sea  level.  Fig.  5b 
gives  the  correspondent  wind  field  which  is  represented  by  iso¬ 
lines  of  constant  velocity.  This  gives  a  more  precise  view  of  the 
wind  field  velocity  in  the  horizontal  plane  (x,  y)  than  the  represen¬ 
tation  of  velocity  given  by  Fig.  5a.  A  maximum  value  of  the  velocity 
reached  6  m/s  in  the  valley  at  x  =  366  km.  It  was  observed  that  wind 
direction  reversal  occurs  when  the  wind  shifts  across  the  valley. 
This  outcome  is  also  predicted  by  Finardi  et  al.  [13]  for  a  relief 
with  a  single  valley.  The  results  obtained  in  our  study  for  differ¬ 
ent  horizontal  planes  show  that  the  rates  decrease  gradually  as  the 
altitude  increases  and  become  regular  at  an  altitude  of  1 000  m.  This 
is  consistent  with  the  results  obtained  by  Sherman  [19]. 


3.3.  Wind  speed  adjusted  values 

The  numerical  technique  was  applied  on  the  model  as  shown 
in  Fig.  3.  The  study  terrain  was  very  rough,  an  elevation 


3.3.2.  Wind  field  velocity  in  vertical  planes 

Fig.  6a-e  show  the  velocity  distribution  in  a  vertical  plane  sec¬ 
tion  of  the  computational  volume  and  in  particular  the  vertical 
motion  resulting  from  terrain  features  and  the  wind  shear  present 
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Fig.  5.  Adjusted  horizontal  wind  field  in  the  ( x ,  y )  plan  above  the  topography. 


in  the  vertical  profile.  The  field  velocity  decreased  when  the  alti¬ 
tude  increase  as  predicted  by  Equation  (12);  the  highest  velocities 
are  located  at  the  wall  (flank)  of  the  mountain  and  the  area  of  the 
valley  located  in  the  intervals  [364,366]  [352,356]  as  described  in 
section  of  horizontal  velocities.  This  area  indicates  the  location  of 
highest  power  wind.  These  results  are  similar  to  those  predicted  by 


the  two  main  diagnostic  models;  MATHEW,  developed  by  Sherman 
[19]  and  MINERVE,  developed  by  Geai  [24]  which  appear  to  have 
been  the  more  accurate  and  the  most  extensively  used. 

Table  2  shows  the  comparison  of  wind  speeds  (m/s)  estimated 
using  the  current  model  and  that  estimated  using  Wasp  for  the 
measuring  station  area  Tenes  [41  ].  The  current  model  gave  a  value 
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Fig.  6.  Adjusted  wind  field  velocities  in  the  vertical  (x,  z )  plane. 


Table  2 

Comparison  of  wind  speeds  estimated  using  the  current  model  with  that  estimated 
by  Wasp  [41  ]. 


Coordinate  of  measurement  point  (m) 

Present  study 
(m/s) 

Wasp 
(m/s)  [41] 

Error 

(X=  3774356,  7=356412.065,0) 

3.47 

3.69 

-6% 

of  3.47  m/s  whereas  the  Wasp  method  gave  a  slightly  higher  value 
of  3.69  m/s  [41  ].  It  is  noted  that  the  result  of  the  estimate  based  on 
the  current  model  is  lower  by  6%  compared  with  that  predicted  by 
Wasp. 


4.  Concluding  remark 

A  numerical  algorithm  for  the  calculation  of  flow  fields  was 
developed.  The  study  showed  that  basic  optimization  of  constraint 
wind  flows  using  a  mass-consistent  model  is  an  effective  method 
for  obtaining  a  3-D  wind  simulation  over  complex  terrain.  An  appli¬ 
cation  was  successfully  adopted  for  a  complex  terrain  in  the  region 
of  Tenes  in  Algeria.  The  current  model  gave  wind  speed  values  that 
were  6%  lower  than  those  estimated  by  Wasp  software.  The  method 
allowed  us  to  determine  the  windiest  areas  in  the  case  study  region. 
The  precision  in  modeling  complex  terrain  and  the  application 
of  the  scheme  difference  method  on  the  topographic  boundary 
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suggest  a  significant  influence  on  the  wind  speed.  The  model  devel¬ 
oped  can  also  be  applied  to  dispersion  of  pollutants,  and  prediction 
of  fire  spreading. 
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